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Abstract 

We propose a numerical variational method for three-dimensional (3D) classical 
lattice models. We construct the variational state as a product of local tensors, and 
improve it by use of the corner transfer matrix renormalization group (CTMRG), 
which is a variant of the density matrix renormalization group (DMRG) applied 
to 2D classical systems. Numerical efficiency of this approximation is investigated 
through trial applications to the 3D Ising model and the 3D 3-state Potts model. 
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1 Introduction 



The density matrix renormalization group (DMRG) [1,2] has been widely ap- 
plied to one-dimensional (ID) quantum systems and two-dimensional (2D) 
classical systems [3,4]. A frontier in DMRG is to extend its numerical algo- 
rithm to higher dimensional systems, chiefly for 2D quantum and 3D classical 
systems. As far as the finite system algorithm is concerned, decomposition of 
higher-dimensional clusters to ID chains proposed by Liang and Pang works 
efficiently [5]. On the other hand, we have not obtained any satisfactory answer 
to extend DMRG toward infinite-size systems in higher dimension. Nishino and 
Okunishi proposed a way of extending DMRG to 3D classical systems, which 
they call 'the corner tensor renormalization group (CTTRG)' [6], as a 3D gen- 
eralization of both the transfer matrix DMRG [3,7] and the corner transfer 
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(a) (b) 

Fig. 1. Two representative constructions of the tensor product variational state: 
(a) the IRF type and (b) the vertex type. The white circles denote spin (or field) 
variables, and the black squares denote auxiliary variables. We consider the simplest 
example of the (a)-type product state in this paper. 

matrix renormalization group (CTMRG) [8,9] for 2D classical systems. Two 
major problems are found in CTTRG when it is applied to the 3D Ising model. 
One is that the calculated transition temperature T c is much higher than one 
of the most reliable T c obtained by the Monte Carlo (MC) simulations [10,11]. 
The other problem is the very slow decay of the the density-matrix eigenvalues 
[12], that spoils the numerical efficiency of the block-spin transformation. 

A way of generalizing DMRG to higher dimensions is to investigate the varia- 
tional structure of DMRG, where the variational state for the transfer matrix 
or the Hamiltonian is written in a product of orthogonal matrices [13]. Two 
types of 2D tensor product states have been considered as higher-dimensional 
extensions of this matrix product state. One is 'the interaction round a face 
(IRF)' type product states shown in figure (1(a)) [14], and the other is the 
vertex type one in figure (1(b)) [15,16]. For a 2D tensor product state V, 
the variational energy and the variational partition function, respectively, is 
written as 

_ (V\H\V) (V\T\V) 

(v\v) and (V\V) ' (1) 



where H and T denotes a Hamiltonian for a 2D quantum system and a transfer 
matrix for a 3D classical system. Let us call such a variational estimation 
as 'the tensor product variational approximation (TPVA)' in the following. 
Calculation of A have been performed by MC simulation [15], product wave 
function renormalization group (PWFRG) [16] or CTMRG [18]. A key point 
in TPVA is to find a good variational state V. So far, they assumed a specific 
form of V, which contains several variational parameters, and tried to find 
out the best V by way of the parameter sweep. For example, Okunishi and 
Nishino [18] investigated TPVA for the 3D Ising model, assuming V in the form 
of the Kramers- Wannier (KW) approximation [17]. Their variational state 
contains two adjustable parameters, and the best V is obtained through a two- 
parameter sweep. Such an intuitive construction of V is, however, not always 
applicable; for example, we don't know how to extend the KW approximation 
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for the 3D Potts models [19]. How can we obtain the best V automatically for 
higher-dimensional systems? We find an answer for 3D classical systems. 



In this paper we propose a self-consistent improvement for the tensor product 
state V by way of CTMRG. We choose the 3D Ising model as an example of the 
3D classical systems, and formulate our self-consistent method in terms of the 
Ising model. In the next section, we introduce the simplest 2D tensor product 
state, and give the formal expression of the variational partition function A 
in eq.(l). We then obtain the self-consistent equation for V in §3, considering 
the variation 8X/8V. In §4 we propose a numerical algorithm to solve the self- 
consistent equation. In §5 we check the numerical efficiency and stability of 
this algorithm when it is applied to the 3D Ising model and the 3-state (q = 3) 
3D Potts model. Conclusions are summarized in §6. 



2 Tensor Product Variational State 



We briefly review the variational formulation of TPVA that was used for the 
KW approximation of the 3D Ising model [18]. Let us consider the 3D Ising 
model on the simple cubic lattice of the size 2N x 2N x oo to X, Y and 
Z directions, respectively, where open (or fixed) boundary conditions are as- 
sumed for both X and Y directions. We are interested in the bulk property 
of this model, and therefore suppose that the system size 2N is sufficiently 
large. Suppose that the neighboring Ising spins a and a' have ferromagnetic 
interaction —Jaa'. The transfer matrix T from a 2N x 2N spin layer 



<7\ 1 



Cl N 



Cl 2N 



&N 1 • • • CjV TV 



&N2N 



(2) 



&2NI ■ ■ ■ &2NN ■ ■ 



&2N2N 



to the next layer [a] is then expressed as a product of local factors 



2JV-1 2N-1 



T[°\°]= n n *«=n*. 



(3) 



i=l j=l ij 



where X- • represents the Boltzmann factor for a local cube 
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Fig. 2. Positions of the spin variables. The plaquett spin {<?ij} in eq.(6) consists 
of 4 neighboring spins, where i' = i + 1 and j' = j + 1, and the cube spin {/U^ } in 
eq.(15) consists of a stack of two plaquett spins {cr^} and {cr^}. 
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parameterized by K = J/k B T. We have used the notation i' — i + 1 and 
f = j + 1. (See figure (2).) We define T[fr|er] so that it is symmetric, because 
the symmetry simplifies the following formulation. 

The variational state in TPVA is a uniform product of local tensors. In this 
paper, we focus on the simplest construction of the tensor product state 



/ 



V[a]=UW, lJ =l[W 



*3 



a i'j a i'j' 



(5) 



where the local tensor does not contains auxiliary variables, which are 
shown by black squares in figure (1(a)) [16,14]; to include the auxiliary vari- 
ables is straightforward, but makes the following equations rather lengthy. 
The tensor product state V[a] is uniform in the sense that is position 
independent. The local tensor has 16 parameters, but not all of them are 
physically independent [20] . For the book keeping, let us use the notation 



KJ = 



a ij a n> 



a i'j a i>j> 



(6) 



for the plaquett spin, and write the local tensor simply as W{a i ^}. In the 
same manner, let us write as X{aij \ a^}. (See figure (2).) 
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Using T[ct|<t] and V[a] thus denned, the variational partition function per layer 
is expressed as 



A = 



E V[a]T[a\a]V[a] £ J] W{a tj } <r y } W{a tj } 
[g]M = [glM g 

' eo>d 2 " " EnK-}) 2 



<7 <7 «J 



where we have defined G° and G 1 as 



(7) 



G°K}=Kki}) 2 , 

G'fal o l3 } = WfajXfal a l3 } W{a i5 } . (8) 

It should be noted that Z° is a partition function of an IRF model [20] on 
2N x 2N square, and Z l is that of a 2-layer lattice model of the same size. 



3 Self-Consistent Relation for the variational state 



Now we explain the self-consistent equation for the variational state V[a], the 
equation which is satisfied when A in eq.(7) is maximized. Let us consider the 
variation of A with respect to the variations of local tensors 



under the condition that the system size 2N is sufficiently large and the bound- 
ary effect is negligible. Then most of the terms in the r.h.s. are almost the 
same, and it is sufficient to consider the variation of A with respect to the 
local change W NN — > W NN + 5 W NN at the center of the system, where W NN 
represents the local tensor at the center. (See eqs.(2) and (5).) 

The variation 5 X/S W NN can be explicitly written down by use of two matrices. 
One is the diagonal matrix 

^{<W = E II (io) 
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where J2[ a y denotes spin configuration sum for all the spins in the layer [a] 
except for the spin plaquett {c NN } at the center; we interpret A{a NN } as 
a 16-dimensional matrix M{a NN \a NN } where M{<t nn \<7 nn } = A{a NN } and 
is zero when {cr NN } ^ Wnn}- F rom the definition, Z° in eq.(7) is equal to 
T,{a NN } G°{cr NN }A{a NN }. The other matrix is 

B{a NN \a NN } = X{a NN \a NN } ^ JJ G 1 {a ij \<J ij } , (11) 

la]'[a]' (ij)?(NN) 



which is related to Z 1 as Z 1 = E{a NN }{a NN } W Wiv) 5 Wwl^JVivl^WTvl- 
By use of A{a NN } and B{a NN \a NN } thus created, we can write down A as 

J2 W{a}B{a\a}W{a} 

x= *em , (12) 

w 



where we have dropped the subscripts from {c NN } and {cxtvat} for book keep- 
ing. The condition 5 \/SW{a NN } = draws the eigenvalue problem 

Y,-r^r B{a\a}W{a} = XW{a} (13) 



{ 



between the matrix A~ X B and W; here we regard W^jc} as a 16-dimensional 
vector. This is the self-consistent equation that an optimized tensor product 
state V[a] should satisfy. 



4 Numerical Algorithm of the Self-Consistent TPVA 



To use the self-consistent relation eq.(13), we have to obtain A{cr} and B{a\a} 
for very large N. Though it is impossible to obtain ^4{<r} and -B{o"|<r} exactly, 
the CTMRG [8,9] enables us to numerically obtain them very accurately. Let 
us introduce a new notation 
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Fig. 3. Positions of the spin variables in eqs.(16). The four white circles denote the 
plaquett spin {a} or {/j,} at the center, and the black squares denote the block-spin 
variables used in CTMRG. 

which groups a pair of adjacent spins and a^. Using fj,^, we can rewrite 
the stack of two plaquett spins {o"^!^-} as 



{Vij} 




(15) 



X{(Jij\ a^} as X{fiij}, and G 1 {a i j\ a^} as G 1 ^^}. (See figure (2)) We drop 
the subscripts from to write it as {//} when its position is apparent. 

The matrices ^4{o"} and B{fi} = B{a\a} can be expressed as a combination 
of the corner transfer matrices (CTMs) and the half-row transfer matrices 
(HRTMs), that appears when we apply CTMRG to both the denominator 
and the numerator of eq.(7) to obtain Z° and Z 1 [18]. Let us write the CTM 
used for the calculation of Z° and Z 1 , respectively, as C°(£cr£') and C 1 
where £, (, and (' are m-state block spin variables. Also let us write HRTM 
as P°(£,aa' £') and P 1 £') 111 the same manner. Note that C°(£o"£') and 
P (£aa'£') are created from and C 1 (C/iC') an d -f >1 (CA t A t 'C) are from 

G 1 -^}. Combining these CTMs and HRTMs, A{cr} and P{/x} are constructed 
as 

M°}= E P°(diW b Q C\i 2 a b Q P\i,a b a c U) C°(& C Q 
Ci-Cs 

^(C 5 /W 6 ) ^(Ce^Cr) ^(Cr^.Cs) ^(Cs^Ci) (16) 

where the positions of the spin variables are shown in figure (3). In principle, we 
can use and thus constructed to solve the self-consistent eq.(13). 
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Fig. 4. Extension of CTM and HRTM [8,9]. The local factor G is created from the 
improved local tensor in eq.(18). 

To make the self-consistent improvement for more efficiently, we em- 

ploy a numerical algorithm that simultaneously performs the extension of the 
system size in CTMRG and the self-consistent improvement by eq.(13). The 
numerical procedures are as follows: 

(a) Create G°{a} and G 1 {fj,} from X{fi} defined in eq.(5) and the initial W^-fo-}: 



G°{a} = (W{a}) 2 , 
G 1 ^} = G x {a\a} = W{a} X{a\ a} W{a} . 



(17) 



The choice of the initial W^{<r} is not so relevant, since it is improved after- 
ward. 

(b) Create the initial C°(£<7£') and the initial P°{£aa'£ r ) from G°{a}, following 
the standard initialization procedure in CTMRG [8,9]. Also create C 1 (C/iC') 
and P 1 (C/U/i'C') from in the same way. 

(c) Obtain the matrices A{a} and B{fx} = B{a\a} using eqs.(16). 

(d) Improve W^{<r} by multiplying A~ X B 



(18) 



and normalize VF new {(T} so that 



E (w nc M) 2 = i 



(19) 



is satisfied. 

(e) Recreate G°{a} and G 1 {fi} by substituting W new {cr} into eqs.(17). 

(f) Extend P° and P 1 to obtain P c ° xt and P^, respectively, by joining the 
recreated G° and G 1 as shown in figure (4); the numerical details are shown 
in ref. [8,9]. Also extend C° and C 1 to obtain C c ° xt and C^ xt . 

(g) Create density matrices from the extended CTMs, and diagonalizing them 
to obtain the RG transformations £ old a -> £ new and CoidJ" ~* dew where £ 
and ( are m- state block spins. Then apply the RG transformations to P® xt , 

-^ext ^cxt an( i Cexf 

(h) Goto (c), and repeat (c)-(g) to improve W^jcr} iteratively, and stop when 
W^a} reaches its fixed point. 
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Fig. 5. Calculated spontaneous magnetization of the 3D Ising model. 

To summarize, we put three additional steps (a), (c) and (d) to the standard 
CTMRG algorithm. 

5 Numerical Results 

Let us check the numerical efficiency and stability of the self-consistent TPVA 
through trial applications to the 3D Ising model and the ferromagnetic q = 3 
Potts model. Figure (5) shows the spontaneous magnetization (a) at the center 
of the 2N x 2N x oo system, where the curve, cross marks, and triangles, 
respectively, represent the result of the MC simulation by Tarpov and Blote 
[10], KW approximation [18], and the self-consistent TPVA. We calculate {a) 
after repeating the iteration (c)-(g) in the last section for N = 10000 times 
at most, keeping m = 10 to m = 20 states for the block spin variables; the 
convergence with respect to m is very fast, where we obtain almost the same 
(a) for the cases m — 10 and 20. The self-consistent improvement by eq.(18) 
is monotonous in the whole parameter range, and no oscillatory instability is 
observed. The calculated transition point K c = 0.2188 is about 1.3% smaller 
than the MC result K^ c = 0.2216544. 

It turns out that the spontaneous magnetization calculated by KW approxi- 
mation, which gives the transition point K^ w = 0.2180, is quite close to the 
result of the self-consistent (SC) TPVA. This means that the intuitive choice 
of the variational state in KW approximation is actually very good, within 
the simplest product state defined in eq.(5). Inclusion of auxiliary variables to 
the tensor product state is necessary for the further improvement of the ten- 
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Fig. 6. Energy per bond (—5 (a, a')) of the ferromagnetic 3D q = 3 Potts Model. 

sor product variational state. Note that the computational time required for 
the KW approximation is several times larger than the self-consistent TPVA, 
because the former finds the partition function extremum via 2-parameter 
sweep. 

Figure (6) shows the energy per bond E = (—5 (a, a')) of the ferromagnetic 
3D q = 3 Potts model, which is calculated by TPVA keeping m up to 15. The 
self-consistent improvement by eq.(18) is again monotonous, and N = 1000 is 
sufficient to get the converged data; we need smaller N for the Potts model 
than Ising model, because the phase transition of the Potts model is first order. 
The calculated energy per bond jumps from E + = —0.5173 to E" = —0.5933 
at the calculated transition point K c = 0.54956, where the calculated free 
energy of the disordered phase coincides with that of the ordered phase. The 
calculated transition point is about 0.18% smaller than one of the most reliable 
MC result K™ c = 0.550565±0.000010 [21]. The latent heat I = 3(£ + -£~) = 
0.22769 is about 41% larger than the MC result I = 0.16160 ± 0.00047 [21]. 



6 Conclusion 

We have proposed a self-consistent TPVA, which gives the optimized tensor 
product state for 3D classical systems, by way of the self-consistent improve- 
ment of the local tensors. Since the method finds out the best variational state 
without using a priori knowledge of the system, the self-consistent TPVA is 
applicable for various 3D models described by short range interactions. 

To generalize the self-consistent TPVA to 2D quantum systems is a next sub- 
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ject that one might consider. This generalization is not trivial, since we have 
used the specific property of 3D classical systems when we obtain the self- 
consistent equation. 
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